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Abstract. We generalize the method of third quantization to a unified exact 
treatment of Rcdfield and Lindblad master equations for open quadratic systems of 
n fermions in terms of diagonalization of An x An matrix. Non-equilibrium thermal 
driving in terms of the Redfield equation is analyzed in detail. We explain how to 
compute all physically relevant quantities, such as non-equilibrium expectation values 
of local observables, various entropies or information measures, or time evolution and 
properties of relaxation. We also discuss how to exactly treat explicitly time dependent 
problems. The general formalism is then applied to study a thermally driven open 
XY spin 1/2 chain. We find that recently proposed non-equilibrium quantum phase 
transition in the open XY chain survives the thermal driving within the Redfield model. 
In particular, the phase of long-range magnetic correlations can be characterized 
by hypersensitivity of the non-equilibrium-steady state to external (bath or bulk) 
parameters. Studying the heat transport we find negative thermal conductance for 
sufficiently strong thermal driving, as well as non-monotonic dependence of the heat 
current on the strength of the bath coupling. 
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1. Introduction 

One of the main challenges of the many-body theory and non-equilibrium statistical 
mechanics is to understand the properties of relaxation of large interacting quantum 
systems. There are two common approaches to this type of problems. One important 
direction is to try to define dynamics in the thermodynamic limit and to investigate its 
properties with rigorous mathematical methods of operator algebras [1, 2, 3]. However, 
in this context explicit results which go beyond existence proofs are quite limited. A 
second approach is to split a large system into a tensor product of a smaller system 
of interest, and the rest (environment), and trying to eliminate all the degrees of 
freedom of the large, macroscopic environment (see e.g. [4, 5]). This approach, although 
involving a series of approximations, is usually more fruitful for explicit calculations and 
quantitative analyses. We may be interested either in relaxation to equilibrium or non- 
equilibrium steady states, depending on the equal or non-equal values of thermodynamic 
potentials assigned to possibly several pieces of environment - which we shall call the 
baths. Such calculations of the quantitative properties of steady states may be very 
useful, for example in the realm of transport theory [6] as may complement the linear 
response calculations and suggest non-linear response or far-from-equilibrium effects. 

However, to date we have had a very few explicit calculations of non-equilibrium 
properties of open many body quantum systems, and mainly they had to focus on small 
systems with a single or a pair of degrees of freedoms (such as spins, or bosons), see for 
example [7, 8] . The reason is that there has been no theoretical techniques to deal with 
open many-body problems except for the Keldysh formalism of non-equilibrium Green's 
functions, which however can easily get too involved for explicit calculations. Recently, 
two new directions have been proposed, both in the direction of numerical simulation 
and theoretical analysis. Namely, in the context of numerical simulations of open many- 
body systems, time-dependent density matrix renormalization group techniques [9] have 
been demonstrated to efficiently simulate relaxation to steady states with the Lindblad 
master equation [10]. On the other hand, it has been shown [11] that the Lindblad 
equation for general quadratic fermionic systems, for example for XY-like quantum spin 
chains which are mappable to quadratic fermionic systems, can be solved explicitly with 
the technique of canonical quantization in the Fock space of operators - third quantization 
for short. 

In this paper we shall show how the third quantization can be generalized to treat 
quadratic systems with arbitrary Markovian master equations , which is not necessarily 
of the Lindblad form. In particular, we shall focus on the Redfield dissipator in 
terms of which we can simulate simple thermal reservoirs, and thermal driving of the 
system under non-equilibrium conditions. After giving a short account on mathematical 
formulation of Markovian master equations and the basic physical assumptions and 
approximations involved in the derivation - in section 2, we shall in section 3 present 
a short but self-contained generalization of the theory [11]. In addition, we shall 
outline the calculation of dynamical correlation functions in Liouvillean dynamics, and 
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formulate an exact treatment of explicitly time-dependent quantum Liouville problems. 
In section 4 we shall apply our technique to treat an open XY spin chain in the non- 
equilibrium Redfield model. We shall outline several intriguing exact numerical results 
on large spin chains. In particular, we show that recently announced quantum phase 
transition in the open XY chain in the local Lindblad bath model, generalizes also 
to non-equilibrium thermal Redfield model with qualitatively identical characteristics. 
The transition is characterized by spontaneous emergence of long range magnetic 
correlations, and hypersensitivity of the steady state to external system's parameters, 
when the transverse magnetic field drops bellow the critical value \h\ < h c = |1 — 7 2 | 
where 7 is the anisotropy parameter. Furthermore, we analyze in some detail the heat 
transport in XY chain, and find regions of negative differential heat conductance for 
strong thermal driving, namely non-monotonic dependence of the heat current on the 
temperature difference between the baths. 



2. Markovian master equations in non-equilibrium quantum physics 

Decomposing the Hilbert space of the universe into a tensor product H = H s <g> Hh of 
the central system Ti s and the bath (or a set of baths) 7ib (environment), one writes the 
total Hamiltonian as 

H = H S ® l h + l s ®H h + \J2 X v® Y ^ (!) 

where X M , are linear operators over 7i s , and linear operators over 7ib- Note that 
X^, can always be chosen to be Hermitian, so this shall be assumed throughout 
this paper. The Markovian quantum master equation for the time evolution of the 
central systems's density matrix p(t) is derived [4] using three main assumptions: (i) 
weak coupling (assuming A to be small), (ii) factorizability of the initial density matrix 
p s (0) <S> Pb(0), and (iii) Born-Markov approximation which rests upon the assumption 
that the bath-correlation functions 

rg^t) := A 2 tr(Y M (t)Y,e-^)/tre-^, %{t) := e^e^ (2) 

decay on much shorter time scale than the central systems dynamics X^(t) : = 
e itH s x^ e -itH s ^y e uge un jt s j n w hich Planck's constant h — 1, and may use different 
inverse temperatures j3 for different pieces of the environment (for different baths). The 
resulting master equation is referred to as the Redfield equation 

^ p (t) = -i[H s ,p(t)]+Vp(t), (3) 
where the dissipator-map has a memoryless kernel with the following general form 

POO 

Vp = J2 drT^(r)[X,(-r)p,X u ] + h.c. (4) 

If one additionally assumes the so-called rotating wave-approximation, one arrives at 
the dynamical semi-group which manifestly preserves the positivity of density matrix 
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at all timesj and can be generally described by the dissipator in the Lindblad form 

V'p = J2^A X ^ X u] + h.c, (5) 

where the only condition is that 7 is a Hermitian 7^ = 7* and positive definite 
matrix. The standard Lindblad form is obtained by diagonalizing the matrix 7 whose 
eigenvectors yield the usual Lindblad operators. The important property of the bath- 
correlation functions (2) (which constitute all that we need to know about the baths) is 
the Kubo-Martin-Schwinger(ifM5) condition 

r^H-i/3) = r^(0, (6) 

which is needed to prove that the thermal state p g ibbs = e~^ Hs / tre" ~" H * is a steady 
state of the master equation (3), provided that all baths are thermalized to the same 
inverse temperature§. However, in case of several thermal baths with possibly different 
temperatures we may expect that p(t) relaxes to a physically very interesting non- 
equilibrium- steady- state (NESS). 



3. Diagonalization of quantum Liouvilleans for quadratic fermi systems 

In this section we give a short account on the general technique of canonical quantization 
in the Liouvile space ('third quantization') and complete diagonalization of Markovian 
master equations (3), with (4) or (5), for quadratic fermionic problems. We treat a finite 
problem with n fermionic degrees of freedom, described by 2n anti-comuting Hermitian 
operators Wj, j = 1,2,..., 2n, {vjj, Wk} = 25j t k, in which the Hamiltonian H may take 
a general quadratic form and the coupling operators may be general linear forms: 

2n 

H s = ^2 w i H i,k w k = w- Hw, (7) 
j,k=i 

2n 

X » = J2 x ^ w i = ^ • ( 8 ) 

3=1 

Thus, In x In matrix H can be chosen to be antisymmetric H T = — H. Throughout 
this paper x = (x±, x 2 , ■ ■ -) T will designate a vector (column) of appropriate scalar 
valued or operator valued symbols Xk- This formalism is immedately applicable either 
for describing, (i) physical fermions c m , m = 1,2, ...,n, where «j 2m -i = c m + cj^, 
ty 2m = i( c m — 4«)> or (") XY-like systems of spins 1/2 with canonical Pauli operators <? m , 
m = 1,2, ... ,n, where the fermionic operators are represented by the famous Jordan- 
Wigner transformation 

W 2m -1 = <7* Yl a m' ' W 2m = <7 y m JJ <C • (9) 

m'<m m'<m 

| This is not the case for equation (3,4) which allows for possible breaking of positivity at initial short 
time interval, the so called sleapage time. 

§ With an additional technical condition of neglecting the Cauchy principal value contribution to the 
time integral 4, see the discussion at the end of subsection 3.2 
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3.1. Fock space of operators 

The fundamental concept for our analysis is a Fock space structure over the 4™ 
dimensional Liouville space of operators /C, which density matrix p(t) is also a member 
of. From here on, we shall adopt Dirac bra-ket notation for the operator space /C which 
is fixed by the following definition of the inner product 

(x\y) = tTx*y, x,yeK. (10) 

We note that 2 2n operator-products \Pa), labelled with a binary multi-index a 

P ai ,a 2 ,..,a 2n := 2~ n ' 2 w?w? • a, G {0, 1} (11) 

constitute a complete orthonormal basis of /C with respect to an inner product. 

In fact it is easy to show that \P a ) is a fermionic Fock basis, and powers 1 in the 
product (11) can be considered like a sort of Fermionic excitations, if we define the 
following set of linear annihilation maps &j over|| IC 

Cj\Pa) = CXjlWjPa), (12) 

and derive the actions of their Hermitian adjoints - the creation linear maps c\ 

c]\P„) = (l-o^K-P,), (13) 

which satisfy canonical anticommutation relations 

{cj,Cjfc} = 0, {cj,c\} = 5 jtk , j,k = 1,2, ...,2n. (14) 

3.2. Bilinear form of the Liouvillean 

The aim is now to show that the generator of the master equation (3) 

C:= -iadiJ + V (15) 

is in general a quadratic form in these adjoint fermionic maps Cj, a-. In order to see that 
clearly, let us define the left and right multiplication maps over /C 

w^\x) := \vjjx), wf\x) := \xwf). (16) 

Inspecting the actions of w^,wf' on the Fock basis |P„) one arrives at the following 
useful identities 

w)=c,+c], (17) 

wf = V(c J -c]) = -(c J -c])V, (18) 

2n 

V := exp(mAf), and Af := c)cj (19) 

i=i 

| We shall use notation where linear maps over the operator space (in physics literature sometimes 
referred to as "super-operators" ) are designated by ". 



where 
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are a parity map, and a number map, respectively, which count the parity and number of 
the adjoint fermionic excitations (number of factors in (11)). Note that V, anticommutes 
with all Cj,Cj, hence the second equality of (18), and V 2 = t. 
The unitary part of the Liouvillean (15) now trivially reads 

-i ad H s = -iw h ■ Hw L + iHw R ■ w R = -4ic f • He . (20) 

The dissipator (4) can be represented as a map over /C as 

2" poo 

P = EE ^ / d ^'^(-r) (r^(r)4 fc + C(r)4 fe ) , (21) 

H,u j,k=l ^° 

where / (t) is a (real-valued) propagator of Heisenberg dynamics in the closed system 

X^t) = x^ • exp(-iad# s t)«; =: f (t) ■ w, (22) 
which - due to (20) - can be explicitly solved for a quadratic Hamiltonian (7), giving 
fft) = exp(4iHt)^, (23) 

and 

£'j, k \ x ) ■= \[wjX,w k }), C" >k \x) := \[w k ,xwj}) (24) 

are fundamental basis dissipators which using (17,18) evaluate to 

£'j,k = tfw* - w\w) = (1 + V){c)c\ - c\c 3 ) + (1 - V){cfc k - c k c]), (25) 
t" hk = w L k wf - wfwf = (i + V){c\c) - clcj) + (1 - V){c kCj - c fc c}). (26) 

It will prove useful if we express the internal dynamics (23) explicitly in terms of 
eigenvalues and eigenvectors of the Hamiltonian matrix H. Since 2n x In matrix is 
anti-symmetric and Hermitian, its real eigenvalues come in pairs e m , —e m ,j = 1, . . . ,n, 
with the corresponding eigenvectors u m1 u* m1 namely Hw m = e m u m and Hw^ = — e m w^ 
since H* = — H. The eigenvectors may and should always be chosen orthonormal (even 
in the case of degeneracies), meaning 

u r u m = 0, «,■«; = S itm . (27) 

Then the spectral decomposition of the Heisenberg dynamics reads 

n 

W = E ( e ~ 4iem< (^ • «J4n + ^ tmt (^ ■ • ( 2 8) 

m=l 

Note that V± = (i ± V)/2 are orthogonal projectors which commute with all 
the terms (20,25,26) that constitute the Liouvillean (15), [V±,C] = 0, and hence the 
dynamics (3) does not mix the operator subspaces /C^ = V±fC composed of even/odd 
number of fermionic operators. Since we are mainly interested in expectation values of 
even observables, such as currents and densities, we shall in the present paper focus on 
the dynamics in the subspace /C + only, and consider the corresponding Lioivillean C\tc+ 

t + = V+CV+. (29) 
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The extension to the odd parity subspace is straightforward. Collecting the results 
(20,21,25,26) it is now obvious that C + is a bilinear form in cj and Cj. For convenience, 
we define 4n Hermitian Majorana maps a r , r — 1, . . . An 

a 2j -i = -j=(6j + c]), a 2j = -j=(cj - cj), (30) 

and express the Liouvillean as 

C + = a-Aa-A l, (31) 

where the 4n x An complex antisymmetrix matrix A, later referred to as a structure 
matrix, and a scalar A , can be expressed as 

A2j-i,2k-i = — 2iHj,k — Mj t k + Mkj, 
A 2j ^ 2k = iM kJ + iM* k , 
A 2j ,2 k -i = - iMj, h - iM* i? 

A 2h2k = -2iHj, k -M* k + M* kJ , (32) 
A = trM+ trM*, 

where M is a 2n x 2n bath-matrix which can be compactly written as 

M:=^^®^, (33) 

V 

_ POO 

^:=E/ <W/„H- ( 34 ) 

Defining the bath-spectral functions f {^(w) := ^ J^dt r^(t)e- iw * for which the KMS 
condition reads 

= ^r^(c), (35) 

and extending the range of integration in (34) to [— oo, oo], or better to say, neglecting 
the Cauchy principal value parts in the integrals - which exactly amounts to neglecting 
the Lamb-shift Hamiltonian term [4] in the master equation - we obtain a very simple 
expression (involving only finite sums) for the bath-vectors 

n 

^ = n J2Yl f ^( 4e ™) ' <^ra + e 4em/3 (^ • ujul) . (36) 

fi m=l 

At this point a remark on neglecting the Lamb-Shift term is in order. As the 
Redfield model already involves a series of physical assumptions and approximations it 
is somewhat difficult to argue under what conditions these terms can be dropped on the 
same level of approximations. However, one can straightforwardly show using the KMS 
condition (6) and Hermiticity (r^ (r))* = (r) that only if the Cauchy principal 
value terms are dropped (i.e. if the range of integration in (4) is extended to [— oo, oo]) 
the Redfield dissipator annihilates the Gibbs state T>\e~ l3Hs ) = 0, and hence Gibbs state 
is the steady state of equilibrium thermal Redfield model. 

Note again that the inverse temperature in (36) could in principle be a function 
of the bath-index (3 = (3 U in case one would be interested in non-equilibrium situation 
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with couplings to several different temperatures. But we should stress that different 
temperatures only make sense among uncorrelated baths for which = for any f3. 

We note also that the present formalism uniformly covers both the Redfield and 
the Lindblad master equations, as the Lindblad dissipator (5) is obtained from (4) by 
simply taking the limit (i) = j^S^t + 0), and then the bath-matrix reduces to a 
Hermitian form M = E^T^®^ = which is equivalent to the one used in [11]. 

3.3. Static Liouvillean: normal modes, non- equilibrium steady state and decay spectrum 

Having the compact form of the Liouvillean (31) - and assuming for the time being that 
the structure matrix A is static i.e. there is no explicit time dependence in the matrix H 
or coupling vectors x^ - we follow Ref.[ll] and explicitly construct its normal form, the 
NESS which is exactly the right-vacuum state of (31) £+|NESS) = 0, the spectral gap, 
and the full spectrum of Liouvillean decay modes, all in terms of spectral decomposition 
of 4n x 4n matrix A. We state the main results here in a compact form. 

Assuming the structure matrix is diagonalizable, its eigenvalues can be paired as 
f3j, —f3j,j — 1, • • • , 2n, assuming Re (5j > 0, and its eigenvectors v 2 j_ 1 (corresponding to 
f3j), and v 2 j (corresponding to — f3f) can always be normalized - irrespective of possible 
degeneracies of among flj, which shall be called rapidities - such that 



VV T = J, J := a x <g> 1 



2n 



/0 1 

10 

1 

10 



\ 

(37) 

/ 



\: ; ; ; 

where V is 4n x 4n matrix whose rth row is given by v r , V TjS := iv, s . Thus the structure 
matrix allows the following decomposition 

A = V T diag{A, -A, . . . , p 2n , -&JJV, (38) 

which after plugging into the Liouvillean (31) immediately brings it to a normal form 

2n 

£+ = -2^&;4 (39) 
i=i 

where 

bj ■= ■ a, b'j ■= v 2 j ■ a, (40) 

are the normal-master-mode (NMM) maps, satisfying almost canonical anti- 
commutation relations 

{b 3 ,b k } = 0, {b J ,b' k } = 5^ k , {b' p b' k } = 0. (41) 

The map bj could be interpreted as an annihilation map and b'j as a creation map of 
jth NMM, but we should note that bj is in general not the Hermitian adjoint of bj. The 
right-vacuum is now essentially defined by 6j|NESS) = 0, whereas the left- vacuum is 
trivial (1|£+ = and satisfies (1|S'- = 0. 
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Assuming that the whole rapidity spectrum is strictly away from the real line 
Re Pj > 0, we state the following exact results: 

(i) |NESS) is unique. 

(ii) Almost any initial density matrix relaxes to NESS with an exponential rate 
A = 2 min Re f3j (the spectral gap of the Liouvillean) . The complete spectrum of 4™ 
eigenvalues of C + is obtained by all possible binary linear combinations Aj, = —2v_-(3, 

^e{0,l}. 

(iii) The expectation value of any quadratic observable WjWk in a (unique) NESS can be 
explicitly computed as 

(wjW k ) NESS := trwjWkp^Ess = 2(l|a 2j _ia 2fe _i|NESS) (42) 

2n 

= 2 ^ f 2m,2j-lf 2m-l,2fc-l (43) 
m=l 
1 f°° 

= --/ du>G 2j -i,2 k -i(w), (44) 

* J-oc 

where 

G(cj) := (A-icul)- 1 (45) 

is a matrix of the non-equilibrium Green's function. The first equality is proven in 
[11] % whereas the last equality requires a simple contour integration on the spectral 
decomposition of the resolvent (45). 

(iv) The Wick theorem may be used for calculation of expectation values of arbitrary 
higher order (even!) observables by sums of all possible pairwise contractions of 
the form (42). 

Note that as soon as some of the rapidities condense to the imaginary axis, or vanish, 
NESS typically becomes non-unique (see Ref. [13] for a detailed discussion of Liouvillean 
degeneracies). 

3.4- Static Liouvillean: time- dependent correlation functions 
The complete Liouvillean propagator can be written explicitly as 

exp(t£ + ) = J2 exp(-2fi/ • P)(b[r ■ ■ • (6 2n )^|NESS)(l|(6 2n )^ • • • {b^. (46) 
^e{o,i} 2 ™ 

It may be of some physical interest to evaluate dynamical response after perturbing 
the NESS by multiplying it with some local observable. In order to avoid discussion 
of negative parity dynamics £_ we take a pair of simplest even-order, quadratic 
observables, and define the corresponding non-equilibrum time-dependent correlation 
function - or non-equlibrium response function - as 

C(j,k),(i,m)(t) ■= (w j (t)wfc(t)w;(0)w m (0)) NESS = 

= 4(l|a 2 j-ia 2 A : _iexp(t£ + )a 2; _ia2 m _i|NESS). (47) 

Small simplification has been made with respect to the statement of Theorem 3 of Ref. [11] which has 
been pointed out by I. Pizorn [12]. 
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Expressing the multiplication maps a 2 j-i = Yll=i(^2r,2j-ib r + V2r-i,2j-iK) an d plugging 
in the propagator (46), while noting that only the terms with or 2 Liouvillean 
excitations contribute, we obtain a simple expression 

(In \ / 2n 

f 2r,2j-lf 2r-l,2fc-l J I ^ ^2r',22-1^2r'-l,2m-l 
r=l / \r'=l 

+ 4 ^ e -2t^ r +M (v 2r ',2j-lV2r,2k-l-V2r,2j-lV2r',2k-l) 
l<r<r'<2n 

X {v2r'-\,2l-1^2r-\,2m-\—^2r-\,2l-\^2r'-\,2m-\)- (48) 



5.5. Time- dependent Liouvilleans 

In this subsecton we indicate how to efficiently treat explicitly time-dependent master 
equations, written in third quantized form as 

^\p(t)) = C + (t)\p(t)), C + (t) = a ■ A(t)a - A (t)l, (49) 

where explicit time-dependece of the structure matrix A(t) may physically arise due 
to driving by means of an external time-dependent force (time dependent matrix 
H(t)) or time dependent coupling operators (time dependent vectors x^{t)). In this 
situation NESS cannot exist, but we shall show that one may still efficiently evaluate 
the propagator 

\p{t))=U\p{G)), U := texp QTdrjC+O-)) , (50) 

where T indicates a time-ordered product. 

The procedure is the following. Note that the space of all anti-symmetric complex 
structure matrices form a Lie algebra so(4n, C). The following straightforward identity 

[|a • Aa, \a • Ba] = |a • [A, B]a, (51) 

holding for any pair of complex 4nx4n matrices A, B, indicates that Liouvilleans (31,49) 
generate 4 n dimensional representation of so(4n, C). Thus, the time-ordered product 
(50) can be evaluated within a Lie group SO(4n, C) of An x An matrices, 

U = t exp ^2 1^ dr A(r)^ (52) 

and 4 " then full Liouvillean propagator is written as 

1 /■* 

U = exp(a-Ca-C l), C = -lnU, C = drA (r). (53) 

2 Jo 

The logarithm of W can be now considered as a 'static' Liouvillean, so we can diagonalize 
it by the methods of subsection (3.3), leading to spectral decomposition of the form (46). 

+ Even if this has to be done numerically, using Trotter-Suzuki decomposition schemes, the 
computational compexity is only polynomial in n. 
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4. XY spin chains 

The theory of the previous two sections shall now be applied to investigate a 
homogeneous, finite XY chain of n spins, described by Pauli matrices <7*' y ' z , j = 1, . . .n 
with the Hamiltonian 

H = J2[ + - V 3 ^ 1 + E ha >> (54) 

which is described by two real parameters, anisotropy 7 and transverse magnetic field /t. 
Without loss of generality we may assume that 7, h e [0, 00). We decide to couple XY 
chain thermally only at its ends, so we consider the most general four coupling operators 
which allow for an explicit solution 

X 1 = Ki((Ji cos#i + of sin^), X 3 = k 3 (<7^cos0 3 + cr^sm0 3 ), 

X 2 = K 2 (a^ cos #2 + °~i sin 6> 2 ), X 4 = K 4 (cr^ cos 6*4 + a Y N sin #4), (55) 

and fully decorrelated baths T@ — S^T^. We take standard baths of harmonic 
oscillators at two ends with possibly different inverse temperatures, and Ohmic spectral 
functions 

rfcM = _ p A,2^/3l, /3 3 ,4^/3r. (56) 

Note that frequency cutoff in the spectral function is irrelevant as we neglect the Lamb 
shift term in the master equation. 

The enitre problem can be fermionized by means of Jordan- Wigner transformation 
(9), namely the Hamiltonian and the coupling operators transform to 

n-l ^ _ ^ n 

if = -i^ (^^"^j^j+i ^-i^-i^S-i^i, 

Xi = cos^i + u? 2 sin^i), X 3 = WK 3 (w 2n cos #3 - w 2n -i sin# 3 ), (57) 

X 2 = /t 2 (u>icos# 2 + u> 2 sin6> 2 ), X 4 = W K A (w 2n cos 9 ^ — w 2n -i sin0 4 ), 

where = (— \) n ~ l WiW 2 ■ ■ ■ w 2n is an operator which commutes with all the elements 
of /C + (or anti-commutes with all the elements of /C _ ) and satisfies WW^ = W^W = 1, 
hence it has no effect on the dissipator (4) in t + . We note however, that the 
commutation of W thru p in (4) for the dynamics in Kr produces a minus sign in 
all the bath terms, i.e. it changes the sign of V-VV-, with respect to a pure fermionic 
problem. 

The 4n x 4n structure matrix has now a specific block-tridiagonal + block-bordered 
form, 

A = A' + B, (58) 
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with 





/a b 


... 


o\ 




cab 


... 





A' = 


Oca 


b ... 







... 


c a 


b 




\0 ... 


c 


a/ 


where a, b, c are 4x4 matrices 






a = —\ht 2 <8 


) a\ b = 


1 

-1 2 
2 


The sequences 


of 4 x 4 matrices \j, l'j 


> r i' 



/ li 


h ••• 




... 




... 


V ri 


r 2 • • • 



ln-1 





r n -i 



\ 



(59) 



(ler — 7(7 



-b J 



(60) 



be straightforwardly computed [seeing (32)] from the form of the coupling vectors 



x 



1,2 



(«i 2 cos 6*1,2, «i,2 sin lf2 , 0, . . . 0)' 



x 



3,4 



(0, . . . , 0, -k 3>4 sin 6> 3>4 , K 3i4 COS 6> 3>4 ) 



and their bath-transformations (36) with (56). Although we are unable to give closed 
form general expressions, we can make an asymptotic estimate - for large n - on the 
decay of these matrices with their distance from the diagonal 



II; 



II' I 



l n+l-j\ 



L n+l-jl 



oc exp(-Kj). 



(61) 

The coefficient K > in general depends only on j,h, and /3l (for LJ or /3r (for r^). 
Note that for the special case of local Lindblad coupling (5) with the same local coupling 
operators (57) , the only non-vanishing blocks which remain are the diagonal ones li 
and r n , given explicitly in Ref.[ll]. 

Below we shall present some intriguing numerical results of the non-equilibrium 
thermal Redfield equation (3,4) for the open XY chain given by (54,55,56), in comparison 
with the local non-equilibrium Lindblad model (5) where a suitable set of coupling 
operators of the form (55) and 4x4 coupling matrix 7 MjV can be chosen to parametrize 

the Lindblad operators Li ;2 = \Jv\ t2 Oi , L 3i4 = ^r^ 2 cr^, parametrized exactly in the 
same way as in Refs.[ll, 14]. For all the numerical results reported for the thermal 
Redfield model we consider the bath parameter values K\ — K3 — 1, k 2 — k± — 0, 
Q\ = 63 = tt/6, and /5l = 0.3, /3r = 5.2 unless /3's are varying, and A = 0.1 unless A is 
varying, whereas for the Lindblad model we always take the bath parameters T\ = 0.5, 

= 0.3, rf = 0.5, = o.i. 



4-1. Non-equilibrium phase transition 

In Ref. [14] an intriguing suggestion of a quantum phase transition far from equilibrium 
in the steady state of an open boundary driven XY spin chain has been put forward. 
Numerical and heuristic theoretical evidence has been given for the spontaneous 
emergence of long range magnetic order in NESS as soon as the magnetic field drops 
below the critical value \h\ < h c , 

/i c = |l- 7 2 |. (62) 
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However, that study was done with local Lindblad reservoirs, so the questions remained 
whether the effect persists in the presence of local thermal reservoirs satisfying KMS 
conditions for non-vanishing temperatures. It is an easy task now to follow the recipes 
of subsection 3.3 and numerically evaluate the spin-spin correlator (note the use of Wick 
theorem as the spin-spin correlator is of fourth order in Wj): 

C />m = tr (of o-^pness) - tr (of p NE ss) tr (o^Pness) (63) 

First, we use efficient prescription (43) to compute correlation matrices at non- 
equilibrium conditions /3l = 0.3 7^ /3r = 5.2 and plot them for two different system sizes 
and five different values of h around h c in figure 1. Results look qualitatively identical 
to those for the Lindblad driving, even for other quantities that were investigated 
numerically in detail in [14]. 



^ = 0.73 h = 0.748 h = 0.75 h = 0.752 h = 0.77 




5 10 15 20 

-iogi D ([ c /,j) 



Figure 1. Correlation matrices C/ m , I horizontal axis (left to right), m vertical axis 
(bottom to top), of the non-equilibrium thermal Redfield model of an open XY chain 
for 7 = 0.5 and different field strength h indicated at the figures (note that h c = 0.75) 
and two diffetrent system sizes n (indicated). Bath parameters are specified in the 
text. 

For example, in figure 2 we plot the phase diagram of the residual correlator 
Cres = Y^im^ 2 |C;,m|/ X^m" 1 '^ 2 1; which also reveals possible criticality in the region 
of a large anisotropy 7 > 1 previously not discussed. We note that size dependence of 
the residual magnetic correlator C rcs shows a very characteristic behaviour: namely 

C res oc exp(— 7]n) with 7] > for \h\ > h c or h = (64) 
C res oc l/n for < \h\ < h c (65) 
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Figure 2. Phase diagram for the non-equilibrium thermal Rcdfield model of an open 
XY chain. We plot the residual correlator C rcs against the bulk parameters 7, h. The 
system size is fixed to n — 100 and bath parameters are specified in the text. 

Thus we shall refer to the regime with < \h\ < h c as long range magnetic correlation 
(LRMC) phase*, the regime with \h\ > h c , or h = 0, as non-LRMC phase, and the 
regime with \h\ = h c as critical. Scaling (64,65) is illustrated in figure 3. Exponential 
decay of the C res (n) in non-LRMC phase (64) is consistent with the exponential decay 
of 2-point correlator with the distance between sites C{r) = Y^j-i= r ^hj/^2j-i=r^ ~ 
exp(— £r), as can be qualitatively noted already in the figure 1. However, we demonstrate 
in figure 4 that the exponents £ could in principle be very different between the Redfield 
and local Lindblad models. Futhermore, as for the Linbdlad model the exponents £ and 77 
[of (64)] appear to be equal, for the Redfield model they don't seem to be simply related. 
Analytical estimation of these exponents present a challenge for future theoretical work. 

However, we note that with the thermal driving with Redfield dissipators, the long- 
range-magnetic order disappears when the temperatures of the baths become equal, 
(h — /3r> an d there we recover, consistently, all the properties of the thermal state [15] 
which are most easily numerically reproduced by the method of Ref.[16] , i.e. fast decay 
of correlations for any h and absence of long-range order. For example, it is interesting 
to note how the residual correlator C res (for large n in the LRMC phase) decreases as 
a function of the difference of inverse temperatures A/3 = /3r — j3^, namely numerics of 
figure 5 suggests clearly that C rcs oc (A/?) 2 . 

Heuristic explanation of this non-equilibrium phase transition is rather straightfor- 
ward [14], however its exact proof and also the quantitative dependence of the decay 
exponent 77(7, h) are still lacking. We note that the transition point h — h c is charac- 

* Note, interestingly, that unlike for the local Lindblad driving[14] the XX line 7 = 0, < \h\ < 1, also 
exhibits long range magnetic correlations for the thermal Redfield driving. 
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Figure 3. Residual correlator C res as a function of the system size n for the LRMC 
phase (7 = 0.5, h = 0.2, left plot) and non-LRMC phase (7 = 0.5, h = 0.9, right plot), 
where we compare the non-equilibrium thermal Rcdficld model (red squares) and the 
non-equilibrum Lindblad model (blue circles) with bath parameters as specified in the 
text. The thin lines indicated the suggested behavior X/n (on the left) and exp(— r/n) 
on the right (with the numerical best fit 77 — 1.192 for the Redfield model and 77 = 0.937 
for the Lindblad model). 



1 r" 




10 20 30 40 



Figure 4. Comparing the decay of the 2-point spin-spin correlator C(r) = 
52i-i=r Qjl X)j-i=r ^ ~ ex P(~£ r ) between the non-equlibrium thermal Redfield 
model (red squares) and non-equilibrium Lindblad model (blue circles) for the same 
values of bulk parameters in the non-LRMC phase (h — 1.05,7 — 0.2, n = 200) and 
bath parameters specified in the text. The thin lines indicate suggested exponential 
decays cx exp(— £r) with the exponents £ = 1.635 (fitting the Redfield model) and 
£ = 0.937 (fitting the Lindblad model). 

terized by a simple property of the XY spin chain quasiparticle dispersion relation 

uj{q) = \J (cos q — h) 2 + 7 2 sin 2 q , (66) 

where 6j = u(2irj/n) would be exactly the (positive) eigenvalues of matrix H if periodic 
boundary conditions would be imposed on the closed system. Namely, in non-LRMC 
phase \h\ > h c there exist only a single pair of trivial stationary points q* = 0,7r, 
whereas in LRMC phase \h\ < h c there exist another pair of nontrivial stationary points 
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Figure 5. Residual correlation C rcs versus the (inverse) temperature drop A/3 between 
the left and the right bath, /3 L = 2 - A/3/2, f3 R = 2 + A/3/2, for the non-equilibrium 
thermal Redfield model of an open XY chain in LRMC phase (7 = 0.5, h = 0.3), system 
size n = 100, and the bath parameters specified in the text. The thin line indicated 
suggested |A/3| 2 behavior. 

±q* 7^ 0, 7T, du/dq\ q=q * = 0, which introduces a new non-trivial length scale 1/q* which 
determines typical sizes of correlated regions in the matrix C^ m (see figure 1). Therefore 
this simple non-equilibrium quasi-particle picture predicts mean-field critical exponent 
l/<f ~ \K - ^T 1/2 as h t h c (confirmed in Ref.[14]). 



< 




10 20 50 100 200 500 1000 

n 



Figure 6. Liouvilllean spectral gap A for the non-equilibrum thermal Redfield model 
of an open XY chain. We plot three different cases with: 7 = 0.5, h — 0.8 > h c 
(non-LRMC phase, light blue circles), 7 = 0.5, h = 0.75 = h c (critical regime, dark 
blue squares), 7 = 0.5, h = 0.3 < h c (LRMC phase, black diamonds), whereas the 
bath parameters are specified in the text. Suggested power law decays n~ 3 and n 
are indicated with thin lines. 

The non-equilibrium quantum phase transition can also be characterized by the 
scaling of the Liouvillean spectral gap A(n), namely in the critical regime one expects a 
qualitative increase in the relaxation time 1/A to NESS. Numerical results (see figrure 
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6) suggest the spectral gap of the Liouvillean remains like in the local Lindblad case 
[11] 

A oc n~ 3 for h ^ h c , A oc n -5 for h = h c , (67) 

although we are at the moment unable to prove this conjecture. Also note slight 
fluctuations of A(n) in the LRMC phase as opposed to a smooth power law in the 
non-LRMC phase. 

Long range correlations for \h\ < h c naturally imply sensitivity of NESS to 
tiny variations in system's parameters. For example, one may expect also that local 
observables in NESS will be then sensitive functions of the bath-driving or even bulk 
parameters, such as the magnetic field h. In figure 7 we plot local magnetization in 
the center of the chain s z = (er?; /0 ) versus the field strength h. Indeed, we notice 

z N ™/ 2/ NESS ta ' 

that for \h\ > h c , s z (h) is a smooth function wheres for \h\ < h c , s z (h) becomes rapidly 
oscillating or better to say, fluctuating, function. Even though the amplitude of these 
oscillations decreases with n, the scale of h on which s z (h) fluctuates decreases with n 
even much faster, so we predict that in the thermodynamic limit n — > oo, in LRMC 
phase the local susceptibility ds z /dh would be ill defined. In summary, LRMC phase 
can be characterized by hypersensitivity of NESS to external parameters. 

0.4 
0.3 

0.1 
0.0 

0.0 0.2 0.4 0.6 0.8 1.0 

h 

Figure 7. Hypersensitivity of NESS to magnetic field strength h. We plot local 
magnetization s z (h) = (< J ^y 2 ) NEg g f° r ^ ne non-equilibrium thermal Redfield model of 
open XY chain with 7 = 0.5 and bath parameters as written in the text. Big blue 
(small red) circles represent data for n = 50 (n = 100), whereas vertical line denotes 
the critical value h — h c . 



4-2. Heat transport and entropies 

An important non-equilibrium physical effect which one can investigate more deeply in 
an open XY chain is the heat transport, which has been recently intensively studied in 
quantum spin chains, see e.g. [17, 18, 19, 20, 21] or [22] for a recent review on the topic. 
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Writing the Hamiltonian (54) in the bulk as a sum H = ^ m H m with a two body 
energy density operator 

rr - 1 +7 1-7 

H m = - 1 — W 2m W2m+l + 1 2 W 2m -lW 2m+ 2 

.h m .hm+l /rc>\ 

~ ly^m-l^m - l^-U) 2m+ lW 2m +2, (68) 

one can derive the local energy current 

Q m = i[H m , H m+1 ] (69) 
= i(l - ^ 2 )(w 2m -iw 2m +3 + w 2m w 2m+4 ) 

- 2i/i(l -7)(w 2 m— lW 2m +l + W2 m +2W2m+4) 

- 2i/i(l + 7)(tt) 2m W 2m +2 + W 2m +lW 2m +3), 

which, by construction, satisfies the continuity equation 

(d/dt){H m ) = (i[H,H m ]) + tr H m Vp(t) = ~{Q m ) + (Q m ^). (70) 

The two terms between the two equality signs above correspond to the unitary and 
dissipative term in the master equation (3). The unitary term has been already 
transformed to a simple expectation value using cyclicity of the trace trx[y,z] = 
tvy[z,x], while the dissipative term can be further shown to vanish in the bulk 
2 < m < n — 2 by excercising the cyclicity of the trace again and transforming the 
integrand of (4) to terms of the form tr X^— t)p[X„, H m ] = 0. The RHS expression of 
eq. (70) then follows from the nearest-neighbour locality of the Hamiltonian. Therefore, 
in NESS the expectation value of the current (Qm) NE ss snou ld be independent of the 
position m. By looking at the dependence of the steady-state current on the system 
size we clearly find ballistic transport, namely (Qm) NESS = 0(n°), irrespectively of 
the temperature differences between the baths and bulk parameters of the model (i.e. 
whether being in the LRMC phase, non-LRMC phase, or critical). However, we find 
very interesting dependence of the heat current on the temperature driving, i.e. on the 
two temperatures of the thermal baths. In figure 8 we plot (Qto)ness versus and /3r 
and find a maximum of the current for intermediate driving, namely when one of the 
temperatures is less than one 1//3l < 1 and the other temperature is about 1//3r ~ 20. 
This is a clear signature of negative differential heat conductance which could perhaps 
be related to similar far-from-equilibroum effects recently observed in spin and charge 
transport [23]. 

This behavior can be nicely characterized by computing the Gibbs entropy of 
NESS. Since NESS is completely characterized by quadratic correlations (wjWfc) NESS 
and the Wick theorem, one can adopt the recipe which has been proposed in Ref. [24] 
for computing block entropies (or entanglement entropies) applied to the entire lattice. 
In fact, taking an arbitrary block region A C {1, . . . , n}, one can compute Von Neumann 
entropy Sa(p) = — t^APA log 2 Pa (in base 2), where pa = tr^p is a reduced density matrix 
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and A denotes the complement of A, as 

S A =J2 + with H ^ := ~ x lo S2 x-(l-x) log 2 x (71) 

and ±iiA,- are the eigenvalues of the 2#(A) x 2#(A) part of the correlation matrix 
Bj t k defined by (wjWfc) NESS =: 5j t k + restricted to Majorana operators Wj,Wk 

corresponding to spins from the block A. The same general procedure has been applied 
to thermal (Gibbs) states in Ref.[16]. When taking the maximal block A — {1, . . . , n} we 
obtain exactly the standard Gibbs entropy of NESS. In figure 9 we plot the Gibbs entropy 
S{i,..., n } as a function of two bath temperatures and show that, quite remarkably, the 
regions of large (maximal) heat current correspond to regions of large (locally maximal) 
Gibbs entropy. This is not unexpected as the product of the heat current and the inverse 
temperature difference A/3 may be understood as the entropy production rate. 

Calculation of Gibbs entropy of NESS provides also a nice way of controlling the 
positivity of NESS as a density matrix, since this is by no means guaranteed by the 
Redfield master equation. Indeed we find that for very small temperatures (large /3's), 
or for very strong bath coupling A, the positivity of NESS might be slightly violated 
(red region in figure 9), namely some of the correlation matrix eigenvalues Vj become 
slightly larger than 1 (but in our numerical experience never by more than 10~ 7 or so). 




I La . . d 
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Figure 8. NESS expectation value of the heat current (QttOness versus two inverse 
temperatures /?l and /3r for the non-equilibrium thermal Redfield model of an open 
XY chain with 7 = 0.5, h = 0.9, sistem size n = 53, and bath parameters given in the 
text. Note that the 'shoulders' of maxima, around /3l ~ 0.05, /3r > 1, and with L and 
R exchanged, could be interpreted as negative differential heat conductance. 

We can use the concept of block entropy of NESS to further characterize the non- 
equilibrium phase transition. For example, we may compute the total (quantum plus 
classical) correlations between two halves of the spin chain in NESS as given by quantum 
mutual information QMI I{n) = S , {i v .. in / 2 } + 5 , { n / 2 +i v .., n } ~ 5{i,...,n}- 
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Figure 9. Gibbs entropy of NESS versus two inverse temperatures /3l and /3r for the 
same parameters as in the figure 8. Note that in the red region (of both large inverse 
temperatures), NESS is no longer a density matrix (at least one of the eigenvalues 
becomes slightly negative) hence the Gibbs entropy is strictly no longer defined there. 

Interestingly, we find (see figure 10) that QMI saturates I(n) = O(n ) in the non- 
LRMC phase (for \h\ > h c ), whereas in LRMC phase (for < \h\ < h c ) QMI becomes 
extensive I(n) = 0(n) indicating a drastic enhancement of correlations in NESS. This 
is again very similar to the behaviour of operator space entanglement entropy (OSEE) 
(analized for the Lindblad model in [14]), so one may extend the relationship between 
QMI and OSEE which has been conjectured for thermal states in Ref.[16] to NESS. 




200 400 600 800 1000 

n 

Figure 10. Another manifestation of the non-equilibrium phase transition: Quantum 
mutual information (QMI) of NESS for non-equlibrium thermal Redfield model of open 
XY spin chain. The bulk parameters are 7 = 0.5 and h = 0.9 > h c = 0.75 (lightest 
blue, saturated curve), h = 0.7, h = 0.5 and h = 0.3 (from lighter to darker blue 
curves). Thin red lines indicated the linear growth of QMI for h < h c . 
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In the context of energy transport it is interesting to look at the energy density 
profiles in NESS. In figure 11 we plot the relative spatial fluctuation of the energy density 
f{m) = \(H rn ) NESS -H\/\H\ where H = YJ^=2 (^)ness is the averaged energy 

density. Quite strikingly, we observe a big variation of f(m) from site to site for LRMC 
phase and very smooth (non-fluctuating) behaviour for the non-LRMC phase which 
is characterized with a bulk-constant f(m) which is exponentially small in n. This 
behaviour can again be considered as a manifestation of hypersensitivity of NESS and 
LRMC. 



50 100 150 200 250 

m 

Figure 11. Another manifestation of non-equilibrium phase transition: positional 
fluctuations in energy density in NESS of non-equilirbium thermal Redfield model of 
open XY chain. We plot the relative fluctuation /(to) = \(H m ) NKSS — H\/\H\ where H 
is the bulk average of energy density (H m }-^ESS- Three curves correspond to 7 = 0.5 
and h — 0.7 < h c (black curve), h = 0.75 = h c (dark blue curve) and h = 0.8 > h c 
(light blue curve), while the system size is n — 253. 

At last we check the dependence of the heat current (Qm) NE g S on the system-bath 
coupling strength A. It was recently reported by Karevski and Platini [25] that the 
spin current J m in the local Lindblad model of an open isotropic XX chain 7 = has 
a non-monotonic dependence on A which can be universally described by a formula 
( J m ) NESS = a'X 2 / {b' + A 4 ) where a', b' are some constants. For the anisotropic XY model 
and general non-equilibrium thermal Redfield driving we are unable to derive an exact 
analytic result, however our numerical simulations suggest a very similar behaviour for 
the heat current 

(Q m ) NESS ««A 2 /(6 + A 4 ), (72) 

where a, b are again some constants which may depend on all system's parameters except 
A. This is particulary interesting as in the anisotropic XY model the spin current is 
not even well defined as there is no corresponding conservation law. This behaviour 
is demonstrated in figure 12 where on may also notice small but detectable deviations 
between numerics and the best fit to (72). We note that the error of the fit does not 
decrease but is roughly constant when we increase the system size n. 
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Figure 12. NESS expectation value of the heat current (Qtti)ness versus the coupling 
strength A for the non-equilibrium thermal Redfield model of an open XY chain with 
h = 0.5 < h c (black curve), h — 0.75 = h c (dark blue curve) and h = 1.0 > h c 
(light blue curve), system size n — 200 and other bath parameters as given in the 
text. Note that the full line gives numerically excellent fit to the Karevski and Platini 
formula [25] (Q m ) = a\ 2 /(b + A 4 ) for best fitted parameters a = 0.040, b = 0.0070, 
a = 0.066, b = 0.0076, and a = 0.088, b = 0.0071 for the three cases of h = 0.5, 0.75, 1.0, 
respectively. 




5. Discussion 

The purpose of the present paper was three-fold. Firstly, we have outlined a general 
method for exact treatment of quadratic many-body Markovian master equations. Our 
formalism, which rests upon treating density operators as elements of a suitable operator 
Fock space (or Liouville-Fock space) is quite flexible and allows for explicit solution 
of static and time-dependent quantum many-body Liouvillean problems, for example 
computation of arbitrary physical obsevables in the non-equilibrium steady state, decay 
rates of approach to the steady state, or even time-evolution of the density matrix 
of externally forced systems described by explicitly time-dependent Liouvilleans, all 
with polynomial computation complexity in number of particles (fermionic degrees of 
freedom) . 

Secondly, we have analyzed in detail the Redfield model of thermal baths within 
our framework. In spite of the fact that the Redfield model does not define a 
proper dynamical semigroup, namely it is not guaranteed to preserve positivity of the 
density operator, we have confirmed that steady states typically correspond to proper 
(positive) density operators. Tiny deviations from positivity have only been observed 
in some test cases for very small temperatures or very large couplings to the baths 
(which anyway violate weak coupling assumption). Furthermore, we have shown that 
coupling the central system with several thermal baths of the Redfield type at different 
temperatures produces physically interesting non-equilibrium steady states, for example 
such states which carry non- vanishing heat current. We wish to stress this physically 
obvious but mathematically delicate point with a particular care, as we have found 
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a qualitatively different result for Lindblad-Davies dissipators which generate proper 
dynamical semigroups and satisfy detailed balance condition with respect to Gibbs 
states [26, 5]. Namely when we constructed a Lindblad-Davies dissipator with respect 
to two baths with two different temperatures coupled to two ends of the system (spin 
chain), we have found that the resulting steady state (fixed point of the Liouvillean 
dynamics) is simply some convex combination of two Gibbs states corresponding to 
the bath temperatues, and as such has always zero heat current and cannot represent 
physical steady state. This implies that the secular approximation (sometimes called the 
rotating wave approximation) which is the one-step from the Redfield to the Lindblad- 
Davies bath model prohibits the emergence of the physical out-of-equilibrium steady 
states with currents, therefore the seemingly harmless rapidly oscillating terms in the 
Redfield dissipator may be absolutely essential for non-equilibrium physics. Thus we 
conjecture that the thermal Redfield model is somehow a minimal mathematical model 
which can describe non-equilibrium thermal driving of a (non-self-thermalizing, e.g. 
integrable) open quantum system. 

Thirdly, we have applied our theory to analyze non-equilibrium quanutm phase 
transition and heat transport in an open XY spin 1/2 chain. We have carefully compared 
numerical results for the non-equilibrium thermal Redfield model and the local Lindblad 
model, which has been discussed before [11, 14]. We have found that the phase diagram 
of the non-equilibrium XY model is insensitive to the theory with which we describe the 
baths, and the differences were only quantitative. In particular we wish to stress that 
thermally driven heat current in the XY chains exhibits non-monotonic dependence 
on the temperature difference which may be interpreted as negative differential heat 
conductance. We believe that our numerical results on non-equilibrium open XY chain 
provide a strong motivation for further analytical work. In particular, we believe that 
the block-tridiagonal plus block-bordered structure of the Liouvillean structure matrix 
(58,59) could be explored in combination with the non-equilibrium Green function 
formula for the observables (44,45) to yield explicit asymptotic results for large n. 

Note added: Formally quite similar approach to non-equilibrium quasi-particles has 
recently been developed independently by Kosov [27]. 

We acknowledge financial support by the Programme Pl-0044, and the Grant Jl- 
2208, of the Slovenian Research Agency (ARRS). 

References 

[1] H. Araki and E. Barouch, J. Stat. Phys. 31, 327 (1983); 

H. Araki, Publ. RIMS Kyoto Univ. 20, 277 (1984). 
[2] D. Rucllc, J. Stat. Phys. 98, 57 (2000). 

[3] V. Jaksic and C.-A. Pillct, J. Stat. Phys. 108, 787 (2002); Commun. Math. Phys. 226, 131 (2002); 

W. Aschbacher, V. Jaksic, Y. Pautrat and C.-A. Pillet, Inroduction to non- equilibrium quantum 

statistical mechanics, in Open Quantum Systems III. Recent Developments Lecture Notes in 

Mathematics, 1882 (2006), 1-66. 
[4] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, 

Oxford 2002). 



Exact solution of Markovian master equations for quadratic fermi systems 



24 



R. Alicki and K. Lcndi, Quantum dynamical semigroups and applications (Springer, 2007). 

H. Wichtcrich, M. J. Hcrich, H. P. Brcuer, J. Gemmer and M. Michel, Phys. Rev. E, 76 031115 
(2007). 

I. Sinaysky, F. Pctruccionc and D. Burgarth, Phys. Rev. A 78, 062301 (2008). 
M. Ban, J. Mod. Opt. 56, 577 (2009). 
S. R. White, Phys. Rev. Lett. 69, 2863 (1992); 

U. Schollwock, Rev. Mod. Phys. 77, 259 (2005); 

G. Vidal, Phys. Rev. Lett. 91, 147902 (2003); ibid. 93, 040502 (2004). 
T. Proscn and M. Znidaric, J. Stat. Mech, P02035 (2009). 
T. Prosen, New J. Phys. 10, 040326 (2008). 

I. Pizorn, "Entanglement and Quantum Critical Phenomena: Operator Spaces and Random Matrix 

Theory", PhD Thesis, University of Ljubljana. 
T. Prosen and I. Pizorn, to be published (2009). 
T. Proscn and I. Pizorn, Phys. Rev. Lett. 101, 105701 (2008). 

E. Barouch, B. M. McCoy, Phys. Rev. A 3, 786 (1971); ibid 2140 (1971). 
M. Znidaric, T. Proscn and I. Pizorn, Phys. Rev. A 78, 022103 (2008). 
K. Saito, S. Takesue and S. Miyashita, Phys. Rev. E 61, 2397 (2000); 

K. Saito, Europhys. Lett. 61, 34 (2003). 

F. Heidrich-Meisner, A. Honecker, D. C. Cabra, and W. Brcnig, Phys. Rev. B 66, 140406 (2002); 

F. Heidrich-Meisner, A. Honecker, W. Brcnig, Phys. Rev. B 71, 184415 (2005). 

C. Mejia-Monasterio, T. Proscn and G. Casati, Europhys. Lett. 72, 520 (2005). 
M. Michel, M. Hartmann, J. Gemmer and G. Mahler, Eur. Phys. J. B. 34, 325 (2003); 

M. Michel, G. Mahler and J. Gemmer, Phys. Rev. Lett. 95, 180602 (2005); 
J. Gemmer, R. Stcinigeweg and M. Michel, Phys. Rev. B 73, 104302 (2006); 
R. Steinigeweg, M. Ogiewa, and J. Gemmer, Europhys. Lett. 87, 10002 (2009). 
L. Arrachea, G. S. Lozano, and A. A. Aligia, Phys. Rev. B 80, 014425 (2009). 
A. Dhar, Adv. Phys. 57, 457 (2008). 

G. Bcncnti, G. Casati, T. Prosen and D. Rossini, Europhys. Lett. 85, 37001 (2009); 

G. Bcncnti, G. Casati, T. Proscn, D. Rossini and M. Znidaric, Phys. Rev. B 80, 035110 (2009). 
G. Vidal, J. I. Latorrc, E. Rico, and A. Kitacv, Phys. Rev. Lett. 90, 227902 (2003); 

J. I. Latorre, E. Rico, and G. Vidal, Quant. Inf. Comp.4, 48 (2004). 

D. Karevski and T. Platini, Phys. Rev. Lett. 102, 207207 (2009). 

E. B. Davies, Comm. Math. Phys. 39, 91 (1974); 
R. Alicki, Rep. Math. Phys. 10, 249 (1976); 

A. Kossakovski, A. Frigcrio, V. Gorini, and M. Verri, Comm. Math. Phys. 57, 97 (1977); 

H. Spohn, Lett. Math. Phys. 2, 33 (1977); 
A. Frigcrio, Lett. Math. Phys. 2, 79 (1977). 

[27] D. S. Kosov, arXiv:0907.1045vl. 



